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ABSTRACT 


The operation of synchronous generators with 
only a direct-axis field winding is severely limited at 
leading power factor. A generator with an additional 
winding on the quadrature axis, provided with suitable 
control, shows considerable advantages over the conven- 
tional generator, 

Thisethesis 2s ‘concerned an, the (farst part 
with a comparison between a conventional wound-rotor and 
a divided-winding-rotor machine, and with the advantage 
gained by the use of the divided-winding rotor machine 
in the system. In the second part the feasibility of 
using a suboptimal control law based on the nonlinear 
differential equations of a divided-winding rotor syn- 
chronous generator is discussed, Numerical results are 
presented for single generator-infinite bus problems 


implementing such control. 
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CHAPTER 1 


INTRODUCTION 


Lee Background 


The annual Canadian electricity demand has 
increased at a rate of five percent per year. There seems 
little doubt that this rate of increase will continue in 
the future, perhaps it may even accelerate. For the post 
war period it was recognized that rapid economic growth, 
necessitates the development of large power systems and 
networks of large transmission capacities. As power sys- 
tems grow,they become interconnected in accordance with 
the development plan for the economic zones in which they 
are located, and afterwards power systems in neighboring 
economic zones will be interconnected. Interconnections 
have many advantages. Some of the advantages are caused 
by diversity (daily, seasonal and annual load diversity), 
staggered installation of facilities, nonsimultaneous 
failures, and scheduled maintenance. Thus, a fault which 
disturbs the stability of the system affects the whole 
system and not only the part in which the fault occurs. 
Without adequate coordination, the interconnected system 
may suffer a catastrophic failure. Experience indicates 
that instability or system collapse, usually develops 
after a period of time as a result of automatic control 


operations. In a transient stability study the action of 
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regulators, governors, relays and other automatic equip- 
ment should be recognized. For the high reliability 
required of a modern power system, stability must be en- 
sured and system control must be automatic. To meet these 
high quality demands the use of automatic control theories 


has been necessary. 


1.2 Power System Stability 


The stability of synchronous machines has been 
o£ concern to the utility industry since its early days, 
because the quality of every system, not only power sys- 
tems, is determined by its stability. To say that a power 
system is stable implies that the energy produced by the 
generators is transmitted to the consumers by keeping 
voltage at the consumers terminals and the frequency in 
the system in certain ranges. Although this capability 
has allowed the treatment of the stability phenomenon as 
one overall problem, still it is useful to continue exam- 
ining the nature of the system behavior under three clas- 
sifications of stability, i.e. steady-state, transient 
and dynamic. 

Steady-state stability 

This type of stability is characterised by the ability 
of a system to maintain gradually attained power transfer 


over the system without loss of stability. Here a dis- 
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tinction could be made between the cases of unregulated 
and regulated systems. The instability of the unregu- 
lated generator is reached at the peak of the steady- 
state power angle curve. The regulated system increases 
the steady-state limit. This stability may be called an 
auuilrerar stavrlicy Dy “wiien We understood an operating 
condition of a synchronous generator, which can only be 
maintained with the aid of automatic control. If the 
automatic control fails the generator goes out of syn- 
chronism., 

Transient stability 

While steady-state stability is the ability of 
the system to operate stably at a given fixed loading and 
transmission system conditions, transient stability refers 
to the ability of the system to maintain stability in the 
presence of a sudden large change in load caused by sys- 
tem switching or by a fault. Generally in a power system 
a transient condition occurs whenever it is undergoing a 
quick change from one state to another, and the variables 
defining the state differ greatly from their values in 
the normal steady state (i.e. angular velocity, angular 
acceleration). The duration of a transient condition is 
short but transient Searcy is of great importance 
for systems. The overall merit of an electrical system 


is almost completely determined by its transient behavior. 
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A system which lacks stability in the transient state has 
no practical value. 

Dynamic stability refers to the ability of the 
system to ensure very rapid damping of small oscillations 
occurring during normal operation (because a succession 
of small disturbances, for example, changes of load, occur 
continually, so that, strictly speaking, the system never 
operates in a steady-state) and to make quite certain 
that self-oscillation, i.e., oscillation with increasing 
amplitude, cannot arise. In effect dynamic stability is 
stability in the sense of small signal control system 


Stability. 


1.3  Divided-Winding-Rotor Power System Stability 


The normal range of stable operation of a gen- 
erator under loaded conditions can be extended by the 
action of a voltage regulator on a direct-axis field 
winding. Growing high-voltage transmission networks and 
increasing use of high-voltage cables produce excess 
reactive power which cannot economically be fully compen- 
sated by parallel reactors, and at a time of low power 
consumption a regulator acting on a direct-axis field 
winding has little effect. An additional winding on the 
quadrature axis, provided with a suitable control can 


increase the stable operating region: that is, cause an 
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improvement of the reactive power-absorption limit at any 
load condition especially under conditions of light sey ee 
Another important aspect of the additional field winding 
is its ability to maintain the air-gap flux at a high 
level, following a fault, and hence to maintain a large 
voltage behind transient reactance. The result of this, 


is high transient-power transfer and good transient- 


stability performance. 


1v4 Stability Control 


Rapid automatic control of generator excitation 
and compensation for line inductive reactance (or in gen- 
eral changes in the network) are two principal methods of 
improving the stability limit of a power system. Both 
methods act to increase power transfer between machines. 

The basic equation describing the dynamic 
(electro-mechanical) performance of a generator in a large 


network is 


2 = —_ 
TP det tkeepdy=be P (leds) 


This is called the swing equation in which Py is the tur- 
bine power imput and Pol is the electrical power output 
of the generator. When a generator is subjected to a 


sudden change in loading there is an imbalance between the 
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mechanical imput and the electrical output. The prime- 
mover cannot act instantly as it has an appreciable time- 
lag. By decreasing the governor time constant and by the 
use of high governor gain the stability limit can be 
increased. Another possibility of control arises from a 
study of the power output of a synchronous generator. 


The power in a steady-state 


Eq V 
e 
or in the transient case 
Eq! Wi 
e 


where Eg and Eq! are the voltages behind the synchronous 


reactance Xa and the transient reactance X respectively. 


d! 
By reducing the transfer reactance between synchronous gen- 
erators (the new transfer reactance is Xa = X, - X.) Say 
by increasing the power transfer, either series or parallel 
Capacitance could theoretically be used, but in practice 
only series capacitors have been installed, owing to the 
Capital cost and technical difficulty of making capacitors 
of suitable rating to operate at line voltages of the 


order of 500 KV. Inserting and removing a series capaci- 


tors in a transmission line can remove the oscillatory 


(Eis), 
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transient due to a fault. This may be done by controllers 
of an elementary type which alter the transmission proper- 
ties in a nearly optimal fashion as a function of the sys- 


16,23 are concerned with 


tem stateS. Some recent papers 
the application of optimal control theory in using series- 
Capacitor switching. 

Another method which also uses changes in the 
network is resistor braking. This device is a bank of 
resistors, located near a large generating plant, con- 
nected in shunt with the three-phase bus through a suit- 
able switch that is normally open. The connection of the 
resistor is initiated by an increase of rotor kinetic 
energy over a suitable threshold value, °° 

The first method of rapid control of generator 
excitation has been comprehensively studied in the tech- 


5,6,10,29 Examining equations (1.2) and 


nieal~diterature. 
(1.3) it can be seen that by controlling the voltage Eq 
behind the synchronous reactance, in steady state condi- 
ELON or by scOnLroulaing the wolcage Eq! behind the transi- 
ent reactance, electrical power output can be controlled. 
Considering not only the simple one-machine infinite bus 
system but also the multi-machine system it is found that 


by raising internal voltages the power that can be trans- 


mitted between any two machines or groups of machines can 
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be increased. The dependance of power on the internal 
woleage or the generator is shownvin Fig. 1, 2,03. IPiq. 
1 is the condition when the field voltage is held constant 


(non-controlled case). 
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For the first controlled case (Fig. 2) the transient vol- 
tage Ea! (which is proportional to field flux linkages 
and for some conditions to the field voltage) is held con- 


stant. The maximum power is appreciably greater. 
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Quick-response excitation, which is able to keep the ter- 
minal voltage constant for all values of load will en- 


sure an almost ideally controlled power limit. 


By changing the internal machine voltage, power output is 
changed. Taking some values of transient voltages and 
considering them as constant for some specific region a 
family of power curves can be drawn as in Fig. 4. 
Comparing the three operating regions of a syn- 
chronous generator from Fig. 1, 2, 3 with the following 
diagram, one more feature of excitation can be noted. 
The power limit is increased beyond 90°, which is a very 
important feature for transient stability, since the 
merit of transient stability is based on the first swing 


(if the system is stable during the first swing it is 


~1e3 oid ict nil ake lal diol | 
-as {fiw heel 36 sausev Lis 02 anptecoo spsdtov fsaim 7 


~ 


jimel Yewed! ballowxsne an weed dial a 
es ae 
a 5 : 


& ee 


el gugivo 1rswoq ,spsticy) estidoem apes eda palpasdo ee 7 
brits aopeslev jaetenss3 to pouiey | omg pnirtsT bopasde 
6 aoipas oiRiogae aoe 703 pasianoc ae merit pataebtange | 
<b spit ni 20 max ad BS, eevw teweg do yLtmel 
-aye 6 to anoipex prisaiago easitt Sd3 priasqmo2 . 
phiwelior ert dtiw £ \s vi .B47 mow? trosstensy- suomomdlp 
«bev0R ad. Keo doktstioxe or stom sno onay ars 


10. 


Stable afterwards). 


5, 5 Ey Bay 


Fig. 4 


As already mentioned in part 1.3 during low 
power consumption, the generator, because of increased 
charging current in the high-voltage transmission net- 
work, often has a leading power factor and may need to 
operate beyond the normal steady-state stability limit 
(which is determined by its capability for absorption of 
reactive power). Regulation on the normal direct-axis 
field winding can not extend the stability range. This 
failure of one field excitation is eliminated by two 
field excitation. The improvement in-.stability can be 


seen from Fig. 5. 
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Figs 5 


Originally the main function of an automatic 
excitation regulator, was to maintain, during normal oper- 
ation, a specified voltage across the system bus-bars 
and supply network. The function of excitation control 
has since been broadened considerably. | Control of the 
excitation under fault condition or immediately after 
fault makes it possible to improve the stability of sta- 
tions operating in parallel, in order to stabilize the 
load, prevent a sudden fall of voltage and ensure satis- 
factory starting of induction motors. All these lead to 
some specific requirements which must be satisfied by an 
automatic excitation regulator. The operating capacity 
of the transmission system should be utilised up to its 


limit. A good margin of transient stability and good 
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damping of hunting should be provided. Switching of the 
number of transmission lines, increases or decreases of 
any kind of load must not deteriorate the stability. This 
can be attained if the exciter has a high sensitivity, 
quick response and a high ceiling voltage (positive and 
negative ceiling have to be implemented, because at some 
Moment* aurapidutall,»ol voltage sis ~required) 5 (To do 

this, the methods of measuring the quantities, which actu- 
ate the automatic control must be simple and reliable. 

The time delay due to inertia in the control device and 

in all components of the excitation system must be small. 


It has been shown 2110727729 


that the stability of a syn- 
chronous machine is improved when the regulators respond 
not only to the variations of the controlling signals 
(voltage, current and angle) but also to the derivatives 
of these variations. Many recently published papers?> "31 
dealing with the application of the optimal control theory 
to power system have proved this. Actually the use of 
velocity and accelerating signals and their combination 


specially through an optimal strategy is most effective 


for system stabilization. 


1.5 Scope of Thesis 


In earlier studies” of the synchronous machine 


with two-axis excitation it was proposed that the machine 
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would normally run with excitation on one field winding only, 
and the control on the second winding used whenever a 
severe disturbance occurs. After steady-state conditions 
have been reached, the quadrature excitation would be 
removed, and normal operation restored. In order to further 
improve the performance of a synchronous machine ,not only 
during the transient period but in all operating condi- 
tions,continuous excitation with independent feedback 
system control is required for each winding, with inde- 
pendent feedback systems. The important feature of these 
two independently controlled windings is that the active 
power P is controlled by the quadrature-axis field which 
is usually called the "torque winding" and the reactive 
power is controlled by the direct-axis field which is 
usually called the "reactive winding". In studies of the 
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d.w.r. machine many references deal with a voltage 


regulator on the direct axis and an angle regulator on 


Bhp oe give a detailed 


the quadrature axis. Some references 
study of d.w.r. machines with various regulators on the 
direct and quadrature axes (voltage regulator, angle 

regulator, proportional voltage regulator, proportional 
angle regulator, derivative angle regulator). The scope 
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machine infinite bus system is used. The optimal stra- 
tegy is formulated through a third order equation (derived 
from the nonlinear one-machine infinite bus model). The 
time optimal strategy gives bang-bang control. Propor- 
tional control as a link between theoretical and practi- 
cal solutions is investigated. After that the combination 
of a proportional signal as a stabilizing signal and simple 


voltage and angle regulators is studied. 
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ST. 
CHAPTER 2 


MATHEMATICAL MODEL 


2.1 Synchronous Generator Equations 


A three-phase non salient-pole machine with 
amortisseur windings is considered. Iron loss is included 
in the amortisseur windings. Saturation is neglected. 
Such a machine has seven windings. On the rotor there 
are two field windings and two amortisseurs windings. On 
the stator there are three stator phase windings. The 
mathematical representation of the synchronous generator 
employs Park's*® equations. The equations which give a 
complete description of the transient and steady-state 
operation of a d.w.r. synchronous generator are summarized 
in the following way: 


Direct axis armature flux linkage 
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Quadrature axis armature flux linkage 
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Direct axis field flux linkage 
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Quadrature axis field flux linkage 


Veg = Ae, nite “F Arka Dea - pata a (2.4) 
Direct axis amortisseur flux linkage 
Vexa ~ *tka t£ka * *kka *ka ~ *aka +a vere, 
Quadrature axis amortisseur flux linkage 
Nees = Se a ea + 7 hee ge - Sees es 2700.) 
Direct axis armature voltage 
AR : W 
e, = — py, - Ri, -—y (25570) 
d Wo d d WG 
Quadrature axis armature voltage 
e = ne pW. - Ri_ + ~~ a (2.8) 
gq SG et gq 5 
Direct axis field voltage 
vie sees eae (2.9) 
ah wie eas fa *fd ; 
Quadrature axis field voltage 
Be seolir eet eee (OF10) 
fq Wo fq EQ mer 
Direct axis amortisseur voltage equation 
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Quadrature axis amortisseur voltage equation 


ghagi a (22) 


The self and mutual reactances used in the above 
equations are expressed as the sum of the mutual inductance 


and the leakage inductance. 


Xs = Xg + Xo] 2 eS) 
Xo = te + Xo (2.14) 
Xecg = Xaqg + Xe] (Ge ES) 
Say = Bogs + Key (22016 ) 
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Mechanical equation of motion 

A synchronous generator is an oscillatory sys- 
tem. A study of this system deals with electro-mechanical 
processes, which are associated with the displacement of 
a heavy rotor and the change of stored magnetic energy. 
The dynamics of the synchronous generator with the speed 


damping term included are represented by the differential 


equation 
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PTC? oa Ege Yea (2.19) 
5 d 
The stored energy which is characterised by the 
term By depends on the machine instantaneous rotor angular 
velocity (or instantaneous frequency) and for a given 


machine is may be expressed as 


(2.20) 


For this study the constant value of Boi i.e., Mes will 
be used (at synchronous speed Wo) Actually transient 
frequency excursions go up to +5% of the synchronous fre- 
quency. Neglecting this gives results which are safe, but 
sometimes pessimistic. 

Electrical torque at air gap 


The term Te on the right hand side of equation 


1 
(2.19) is the electromagnetic torque which has to match 
the mechanical turbine torque Ty, and thus prevent torque 


differentials that pull the generator out of synchronism. 


The electromagnetic torque is expressed as: 


Tol = va = - Vag ig (25:2 1)) 
The first component is due to the interaction between the 


direct-axis flux and the quadrature-axis current, and the 


second component is due to the interaction between the 
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quadrature-axis flux and the direct-axis current. 


2026 Single-Machine Infinite Bus*system 


The present study is concerned with the system 
of Fig. 2.1 in which a generator, connected to an infin- 
ite bus through a two circuit line is equipped with a voltage 
regulator on one field winding (direct axis) and an angle 


regulator on the other (quadrature axis). 


Mae = : 


angle 
regulat. 
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To complete the mathematical model of this sys- 
tem equations of the transmission line and the regulators 
have to be added to the equations of the synchronous gen- 
erator. 

Transmission system 

The transmission system is represented by lumped 
series inductance and resistance in Park's reference frame. 
The infinite bus system is characterised by constant fre- 


quency and voltage and zero impedance. 
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Angle regulator 
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where u(t) is a bang-bang or a proportional stabilizing 


Signal. 


ome Steady-State Phasor Representation 


Phasor diagrams are yvery useful in the theory of 
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electric machines. They give relations between the vari- 
ables as the operating conditions change. But the distinc- 
tion between the phasor diagram and stability is an impor- 
tant one. The phasor diagram tells nothing about stabil- 
ity. To get a better picture of the d.w.r. generator the 


conventional phasor diagram for a generator is drawn. 


{q.a. 


Fig. 2.2 is a conventional phasor diagram (excitation on 
direct axis only) at a lagging power factor. Fig. 2.3 
shows the diagram for the same operational condition when 
the machine has a quadrature winding control which can 


hold 6 at a specific value. 
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It was mentioned in part 1.5 that one of the advantages 
of the d.w.r. generator is that the active power P is con- 
trolled by the quadrature axis field. It will be confirmed 
in the following. 

From Fig. 2.3 it can be seen that ve is com- 


posed of two components Val and Voq' 


V..siné - X_I 


Vod B ad 


V.cOss. axa t (2526) 
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The excitation voltage components Voa and Vea (behind syn- 
chronous impedances Xq and a respectively) are directly 


PLOpOLtLonal to'the: currents ote: rleg respectively. 


Viod cs Xagl tq 


ee = X qt ea (Qet2y)) 
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By resolving the stator current into two components i and 


I. the power of the synchronous generator is 


Oo = Wed (2528) 


where I_ is the component of stator current in the phase 


with bus voltage V., and the I. component of stator cur- 


B 
rent lags the bus voltage by 90° (Fig. 2.4). 


Fig. 2.4 


The relation between two sets of current com- 


ponents Torly and tartg are: 
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p gq d 


Hi 
il 


= = i .29 
I. Ia cosé ae sind (2,29) 


“— a a Als 


aia | te , Ws 
i) 7 LM 


(See) 


metic) ‘aovste 2 3p dam yet 


HBS ay ‘0g em 


=a al 7 


24. 


From the equations 2.23 and 2.22 expressions for Se and 


Ig are: 
x V 
Hs = eae a + = sing 
g g 
og V 
nS ae Teg - yo cosé (2.30) 
d d 


Substituting 2.30 into 2.29 and then into 2.28 expressions 


for power of génerator are: 
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Tf the angle 6 us held constant at zero value 
the expressions for the power of a/synchronous generator 


are: 
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The active power ae is controlled only by the quadrature 
axis field winding, while the reactive power Q, is con- 


trolled by the direct axis field current. 


2.4 The State Space Equations of the Single-Generator 


Infinite Bus Case 
Combining equations 2.1-2.12 and 2.22-2.23 gives 


the form 
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tae Spy oes kq 
ele Once lia pac 
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The A and E coefficients in 2.34 and 2.35 are given in 
Appendix I. 

By separating derivatives in equations 2.34 and 
2.35 and defining new coefficients (C, D whose values are 
given in Appendix I), the differential equations of the 
synchronous generator-infinite bus case, in matrix form 


are given as: 
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Ph pie Operating Point Calculation 


An operating point is determined by Por QO wand 


Oo 


Sos The values of Vegr Veg? leg and tq 


and can be calculated using matrix form 2.36. 
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CHAPTER 3 


CLOSED LOOP SUB-OPTIMAL CONTROL OF A D.W.R. 


SYNCHRONOUS GENERATOR 


3.1 General 

Modern optimal control theory has been increas- 
ingly applied, because a growing concern for the economic 
operation of complex industrial and manufacturing processes 
requires a new way of planning and design. Because of 
computer availability, optimal theory as an area of math- 
ematics has made a significant impact on the design and 
Operation of systems, small and large. The usefulness of 
optimization theory is dependent upon the ability to 
obtain numerical solutions. Having numerical solutions 
we can test the validity of conventional systems by com- 
parison. Such a comparison must consider other factors 
such as the availability and cost of components for con- 
structing the optimal control plus the cost and time 


required by the digital computer. 


3-2) etatement sof the Problem 

For the fastest transient removal in the given system 
determine the controls Vea and v.85 which force any given 
initial state (Cys Gor G3) to the stable equilibrium 
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delays of the regulators on the direct and quadrature axes 
are neglected. The magnitude of the field voltages are 
bounded and are piecewise constant function of time (bang- 
bang control). Their signs are determined by the time 
optimal scheme. Now, the new control policy can be inter- 
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The problem can be restated as: 
Determine the control Vea and Veg? subject to constraints 
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3.3 Optimal Strategy 


The swing equation of synchronous generator- 


infinite bus case is given by 
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Differentiating the above expression yields 
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The turbine torque is considered constant in this analy- 
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where all the above derivatives are defined in section 
es tlve 

Substituting expressions for armature fluxes and 
eliminating a few terms because we are dealing with a 


nonsalient generator (i.e., Xq = =e equation 3.4 becomes 
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Further elimination of derivatives in expression 3.5 and 
grouping of terms is given in Appendix II. Using the 


results from Appendix II expression 3.3 is 


3 2 ran 
T5p Sige EG isla alesse T3)\et ie BK) (3.6) 
or 
K Ww 
3 Z6tded2 i. O° 
Dac. ot pth aa tan 0 ae Veg BN + Veg BK) Ga ea) 
Further simplification gives 
Beh a Bes Sa (3.8) 
%6 onl 
where K = z is a positive constant and T behaves as a 
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Defining the variables Y,(t), Yo(t), Y(t) 
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Y,(t) = pd 
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The state space equations of 3.8 are 
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or in a vector-matrix form 
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VACY 2 AVYAC) +°7T (71) 
The matrix A can be transferred to its Jordan canonical 
form because two of three eigenvalues are the same. The 
transformation matrix P is given such the relation J(A) = 


p-1 ap is satisfied, i.e. 
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ihe vector Y(t) can be. transformed into vector Z(t) directly 
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Differentiating this equation yields 


p Z(t) = P+ py(t) (3.14) 


and then eliminating vector Y(t) 
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Bay Sak Z,(t) + T(t) (307) 


For simplicity, new state variables x,(t), X5(t), X(t) 


can be defined by setting 


5) 
x, (t) = K Z, (t) 
Pe 
Xx (t) = K Z(t) 
x, (t) = K Z(t) (Seo) 


After this linear transformation the system of differential 


equations in the variable x are 
px, (t) = K x, (t) 5S MANGE 
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Now, having system 3.19 our problem can be restated. For 

a given system 3.19 find thejadmissible “control that forces 
the system of Fag. 3.1 from any initial state (Si, Gor G3) 
to the stable operating point in minimum time. The gen- 
eral expression for the Hamiltonian H for the minimum-time 


control of the system 3.19 is 
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Equation 3.20 applied to system 3.19 gives 


H= 1+ kK x5 hy ke, h, + K DCE Cer (t) + xe) + A, (t)) 


(377219 


Tne eheMinimal control, 1.e., the controlswhich minimizes 


the Hamiltonian, is given by 
Sanit )i -sgn{-A, (t) + hy (t) + ha (t)} (3,22) 
where the control variable is constrained by the inequality 


min max 


The costate variables Ay (t), hy (t) and ha (t) satisfy the 


equations 
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is a linear combination of constants, linear function and 
one exponential function. This kind of function can have 


at most two zeros as shown in Fig. 3.2. 
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tion 3.19 can be solved using T(t) = const, with the ini- 
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These two equations determine the trajectory in the three- 


dimensional state space (with the initial state oie Soy ba). 
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is the equation of the trajectories which are parabolic. 


The family of these curves is shown in Fig. 3.3 
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For some cases our objective is to drive any initial state 
to the origin of the state plane. The parts of the above 
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the system to the y switch curve (use of the other control 
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where T = ane ante x5>0 and T = ear sh x5<0. The second 


equation of the system (3.26) 
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gives the trajectories in the Koy X3 plane with various 
initial states (Sor G3). ETO, eo". Sle 


Let us consider the two trajectories which pass 
through the oniganlingthe X51 X3 plane. By separating the 


variables and the initial coordinates in 3.28 the result 
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of the X51 X3 plane (i.e. xX, = 0, X3 = 0) two of three 
conditions for equilibrium are satisfied, because xe 0 
and X32 = 0 implies pdé = 0 and p-6 = 0. The third one 
6 < =) is already satisfied. For cases in which 
the rotor pe is larger than Zu, the control scheme 


has to be decided from the combination of the X41 Xo and 
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X, planes, i.e., from the switching function in space 
X11 Xo X3- The projection of the switch curve y in the 


X5 plane is given by 3.26. 
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iki, > eye ee, ES (3.32) 


Let us now consider an initial state (Xi, Xor X3) 
which belongs to the trajectory in the space which inter- 
sects the y curve, and such that a system staying on that 
trajectory can be forced to reach the y curve Fig. 3.6. 
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Eliminate X191 Xo X39 from the set of equations 3.34- 


S.3/ by first finding X55 from equations 3.34 and 3.37 
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Eliminate x from equations 3.35 and 3.37 
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Substituting equation 3.37 into equation 3.39 gives 
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In equation 3.41 there is an ambiguity of sign. 
Going back to the equations 3.36 and 3.37, i.e., to the 
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The control T = Ty eee 2b tmax has to be decided such 
that the trajectory generated by one of them intersects 
the switching surface (Fig. 3.8). For a state above the 
switching surface the control is T = Tee’ and for a state 
below the switching surface the control is T=T.. ., 

In the control scheme we first have to find the 
value of the angle 6. If the angle is out of the range 


- = the switching surface has to be used to bring the 
fe) 


angle 6 into the range +t = -) “Abter, that tne control 
fe) 
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scheme from the plane x x, is decided. 


Deli 
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47. 
CHAPTER 4 


4.1 Method of Simulation 
The equations were solved using a Runge-Kutta-Gill 
4th order integration routine with a step size of 0.001 sec. 
on an IBM-360 model 67 (Appendix III). 
Several forms of disturbance were assumed for 
computation of transient characteristics. 
Disturbance considered: 
a) Torque steps of 10, 30 and 50% of the system load for 
a Single machine infinite bus system. Fig. 2.4 and 
Appendix I. 
b) Torque pulse of 100% for the same system and the same 


operating points. 


A 2X Cita Lon .COntrols 

1. Bang-Bang Control 

From the optimal trajectories and the switching 
function note that the control function T has two values. 
a and Th . The optimal scheme requires that the value 


max in 
of T changes instantaneously from one to the other. T is 


given as 


UW) 
ae ae 
T= ---[E + Veg BN + Ve, BK] 


fq 
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To change the value of T to either Tha or to the i the 


sige) ax 


actual control : ; 
controls have to change from Veg an? Veg ae ee) 
Ved max’ Veg max * 


The field voltages are constrained in magnitude. 
Usually this, constraint is t5/funiets,. 

The angle time characteristics (Fig. 4.1) show 
responses of controlled and uncontrolled machine for a 10% 
torque step. For the uncontrolled case the system stabil- 
izes but is poorly damped. For the machine (Fig. 4.1-b) 
with bang-bang control, the time required to reach a 
stable point is short and the response is heavily damped. 
The phase plane plot (Fig. 4.2) shows the nature of the 
forced and free damped oscillatory systems. For the 
forced system (Fig. 4.2) the rotor velocity is brought 
near a stable equilibrium point in a fairly short period 
of time but exhibits forced oscillations characterised by 
slow variations of speed due to the large rotor inertia. 

The responses for more severe disturbances of 
30% and 50% torque steps are shown in Fig. 4.3 and Fig. 
4.4. As can be expected the final angle becomes bigger and 
the time required to reach a stable point is increased. 
Angle time characteristics from Fig. 4.3 and Fig. 4.4 are 
slightly oscillatory particularly in the initial period. 
To analyse this we have to go back to the switching trajectory. 


Starting from any initial state we want the forced trajec- 
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Figg 4.1 
a) unregulated machine 


b) regulated machine (bang-bang control) 
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Velocity vs acceleration plot corresponding to Fig. 4.1 
a) unregulated machine 10% torque step 


b) regulated machine (bang-bang control) 
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Angle time characteristics for 30% torque step 
a) unregulated machine 


b) regulated machine 
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Angle time characteristics for 50% torque step 
a) unregulated machine 


b) regulated machine (bang-bang control) 
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tory to reach the switching trajectory, and after the 
transition from one control to another occurs, the forced 
trajectory should continue along the switching trajectory 
torthe, origin. 
Due to system imperfections the forced trajectory paral- 
lels the switch curve. For a big disturbance and strong 
control the forced trajectory tends to recross the switch 
curve after having switched once and control action steps 
back and forth causing the trajectory to zigzag along the 
Switching curve. 

A more important cause of oscillatory responses 
is /due to the/"control function. 


Ww 
eae. 


T(t) = 7 


BNICL) + WV 


JE (t) + Veg 


a BK (t)| 


1S a function of the terms E(G),)BN(t)i, (BK(ty (i.e.,ja 
function of state) and the field voltages Veg and Veg: 
From this expression it is seen that the system can be 


controlled only if the relation 


E(t) i (Veg BN (t)) +1 V BK (t) ) 


fq 


is satisfied. By checking the variation of E(t), BN(t) 
and BK(t) it was noted that there were uncontrollable ini- 
tial periods of short duration. 

Fig. 4.5 shows the comparison between the varia- 


tion of field voltages and field currents. The stepwise 
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field voltage variation cannot be followed by field cur- 
rents, because of field time constants. In contrast to a 
Single field machine the field currents may reach negative 
values to produce increased synchronizing (or resynchron- 
izing) torque. 

2. Proportional Control 

Instead of using bang-bang control, i.e., the 
values of field voltages decided by the sign of switching 


function, a proportional control may be used. 


(4.2) 
With constrains on Veg and Ue as defined in 3.2 


The values of K should be as large as possible for K = 


1,2 tye 


©, equations 4.2 reduce to a bang-bang control. 

Unfortunately large Ki, K, may throw the system 
into a steady hunt. 

The comparison between a bang-bang and propor- 
tional control for a 100% torque pulse of 3 cycles duration 
is given in Fig. 4.6. The angle time characteristics in 
Fig. 4.6 show that the bang-bang control has smaller first 
and second overshoots. The final rotor angle is also 
smaller than for the proportional control. 

From Fig. 4.11 it can be seen that the bang-bang 


and proportional control deteriorate the terminal voltage, 
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Fig. 4.6 
Angle time characteristics for 100% torque pulse (3 cycles) 
a) unregulated machine 
b) regulated machine (bang-bang control) 


c) regulated machine (proportional control) 
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Velocity vs acceleration plot corresponding to Fig. 4.6 


100% torque pulse (case b, bang-bang control) 
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ACCELERATION 


VELOCITY 


Velocity vs acceleration corresponding to Fig. 4.6 


100% torque pulse (case c, proportional control) 
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Field voltages and field currents time characteristics 


corresponding to Fig. 4.6 100% torque pulse (case b, bang- 
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because the time optimal policy was defined to bring sys- 
tem to a stable angle in a minimum time. 

3. The Combination of the Voltage Regulator and 
the Angle Regulator with the Proportional Stabilizing 
Signal 

Practical considerations require a constant angle 
value and voltage within a given range. This requires 
implementation of a voltage regulator on the direct axis 
and an angle regulator on the quadrature axis. By intro- 
ducing a proportional stabilizing signal on both axes of 
the rotor the final rotor angle is reached in shorter time. 

Fig. 4,.12-b is an angle time characteristics for 
a 100% torque puls. The new final a teoie angle is reached 
in a short period of time with only one overshoot and one 
undershoot. From Fig. 4.14-b it is seen that the varia- 
tion of terminal voltage is much improved (compare with 
Big.4. Lib bang -bang sCOnNcrol) 

For the same system and for the same disturbance 
an angle regulator on the quadrature axis is added to bring 
the system back to its initial roton, position. Fid. 4712— 
c shows that after the disturbance is cleared the initial 
rotor angle is regained. 

From Fig. 4.14 it can be seen that the terminal 


voltage is initially deteriorated by the strong stabilizing 


Signal. 
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At the initial time,variations of field voltage 
are almost bang-bang (Fig. 4.15) but after the stable point 
is reached the voltage variation on direct axis is smoother 
and on the quadrature axis is still big due to the propor- 


tional signal presented on quadrature axis. 


ANGLE (deg.) 


+20 


Fige, 4 $L2 
Angle time characteristics for 100% torque pulse (3 cycles) 
a) unregulated machine 
b) regulated machine with a proportional control on quadra- 
ture axis and voltage regulator on direct axis. 
c) regulated machine with an angle regulator on quadrature 
axis, voltage regulator on direct axis and stabilizing 


signal on both axis. 
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Figs 4.13 


Velocity vs acceleration plots corresponding to Fig. 4.12 

a) regulated machine with voltage regulator on direct axis 
and proportional control on quadrature axis 

b) regulated machine with voltage and angle regulators 


plus stabilizing signal 
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Fig. 4.14 
Terminal voltages plot corresponding to Big. 4. h2 


unregulated machine 


regulated machine with voltage regulator on direct axis 


and proportional control on quadrature axis. 
regulated machine with voltage and angle regulators 


plus stabilizing signal on both axis 
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Field voltages variation for 100% pulse (3 cycles) corres- 
ponding to Fig. 4.12 (case b, voltage regulator on direct 


axis and proportional signal on quadrature axis. 
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Fig. 4.16 


Field voltages variation for 100% torque pulse (3 cycles) 
corresponding to Fig. 3.12 (case c, voltage regulator on 
direct axis and angle regulator on quadrature axis with 
proportional signal on both axis). After reaching stable 


angle the proportional signal is switched off. 


80 


(awe) ANT 


ent Pope. Jb) | ba: a 
(092) aMiT 


a1. pba 


(eoloyo £) selag s etipxos goo 363, | peadstzsy aepssiov plein 
10 rodsiupe: soraion ‘e 9289) st +Bhe oF paibnoqsexzes, 

ijcw eixs ewigeabeup iat) sotplopex ofpns, bas atxs soexib 
eldate paidoset erin . (abe died no Lanpte Lanoiszogorg— 


bo ila ek ee ieisinei ait elpas 


CHAPTER 5 


CONCLUSIONS AND SUGGESTIONS FOR FURTHER RESEARCH 


opal Summary and Conclusions 


This thesis presents a study of the transient 
stability of a dual-exciter synchronous generator. Opti- 
mal control theory in the closed form was applied to a 
nonlinear model. It can be noted that the unknown initial 
values of the costate variables Tir To and 3 (in equation 
3.22) were not calculated. The sub-optimal control was 
found solely as a function of the states. Such an imple- 
mentation requires the measurement of the entire state 
vector. The result of this is a bang-bang control law, 
i.e., control actions switch between extremes of allowed 
values. Bang-bang control lacks the ability to cope with 
feedback control required for other purposes. As a prac- 
tical solution a proportional control in combination with 
a voltage regulator on the direct axis and an angle regu- 
lator on the quadrature axis is proposed. 

The significance of the analysis of this non- 
linear model is of importance because it gives a correct 
representation of the system in any operating condition. 
Control of such a system is difficult especially with 
single control. The analysis of the results shows that 


the time optimal responses were not achieved. In general 
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sub-optimal solutions were obtained. The time required 
for stabilization was less than one second compared to a 
time of three to four seconds obtained using conventional 
regulators with first and second derivative feedback sta- 
bilizing signals. 

The double winding rotor machine while somewhat 
more expensive than a conventional machine has significant 
operating advantage. Use of sub-optimal closed form exci- 
tation as described in this paper can provide stabiliza- 
tion at a cost which is small compared to the costs of 
presently used alternatives. 

With judicious choice of approximate contra Sig- 
nals the control complexity can be reduced to little more 


than that of conventional excitation control systems. 


5.2 Suggestion for Further “Research 

One factor influencing the performance of the 
system is the uncontrollable initial periods. This depends 
on the severity of the disturbance. Further study on modi- 
fications of the control may reduce or eliminate this. 
This will increase the complexity of the system. Another 
area for further study is that of obtaining simple satis- 
factory practical control signals. That is, a control 


‘using readily measured variables. 
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NOTATION 


Direct-axis field voltage 

Quadrature-axis field voltage 

Direct-axis bus voltage 

Quadrature-axis bus voltage 

Direct-axis armature voltage 

Quadrature-axis armature voltage 

Direct-axis field resistance 

quadrature-axis field resistance 

Direct-axis amortisseur resistance 

Quadr ature-axis amortisseur resistance 

Armature resistance 

Direct-axis armature reactance 

Quadrature-axis armature reactance 

Direct-axis mutual reactance between armature and 
field 

Quadrature-axis mutual reactance between armature 
and field 

Direct-axis mutual reactance between amortisseur 
and field 

Quadr ature-axis mutual reactance between amortisseur 
and field 

Direct-axis field reactance 


Quadrature-axis field reactance 
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76. 
Xekd Direct-axis amortisseur reactance 


*kkq Quadrature-axis amortisseur reactance 


x, Line reactance 

va Direct-axis armature flux linkage 

Yq Quadrature-axis armature flux linkage 

beg Direct-axis field,flux linkage 

Veg Quadrature-axis field flux linkage 

beg Direct-axis amortisseur flux linkage 

Veg Quadrature-axis amortisseur flux linkage 

Wo Synchronous speed in radians per second 
Differential operator 

ce Inertia constant in seconds 

6 Rotor angle with reference to infinite bus in radians 

T, Turbine torque 

Tol Electromagnetic torque 

Ka Mechanical damping coefficient 

Ie Voltage regulator time constant 

T Angle regulator time constant 

K Voltage regulator gain 

K. Angle regulator gain 

deg Direct-axis field current 

Leg Quadrature-axis field current 

ig Direct-axis armature current 

fe Quadrature-axis armature current 


u(t) Stabilizing signal 


3 Switching function 
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The A(I,J) coefficients are: 
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A 2) 
A(1,3) 
A(1,4) 
A(1,5) 


A(1,6) 
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A(4,2) 
A(4,3) 
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A(4,5) 
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A373) 
A(3,4) 
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The C and D coefficients are: 
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Coe =, (Ao Aas ay 223) 72N 
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Cao = (Az) Ayo ~ Air 432)/4u 
Co = (Aqy Ago ~ Aga Ar2)/Au 
Diy = (E52 B33 ~ #32 ®23)/Fu 
Dip = (E32 Fy3 ~ Bio ®33)/Eu 
Dy3 = (Ey2 Fo3 ~ Fo2 13)/Fu 
Doz = (B31 2o3 ~ Bai ®33)/Em 
Do> = (By) 233 ~ #31 #13) Eu 
Do = (Ep) Fy3 ~ ®yi ®23)/Pu 
D3, = (Eg, 232 ~ E31 B22) /Pm 
D35 = (23, Ey. ~ Bui F32)/Py 
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System Parameters 
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16.5 KV turbo- 


generator which has a quadrature-axis field winding iden- 


tical with the direct-axis field winding. 
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The operating point calculated from 2.35 is given by: 
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Substituting expression for derivatives in 
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Equation A2-1 can be written as 
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APPENDIX IIT 


Computer Program 


READ THE INITIALS 
VALUES & PARAMETERS 


COMPUTE COEFFICIENTS OF 


THE MATRIX Y=AY+BT 


CALL DRKGS 
CALL FCT 


COMPUTE E,BN,BK,T1, 
12, Ke XD 1X3 


oo 
a 


COMPUTE THE VALUE OF 
SWITCHING FUNCTION Z 


VFD=+5SIGN(BN) | VFD=-5SIGN (BN) 

VFQ=+5SIGN (BK) VFQ=-5SIGN(BK) 
: COMPUTE Y 
CALL OUTP 


WRITE X,(Y(1I),I=1,8)) 
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FORTRAN IV COMPILER (20.1) MAIN 08-08-72 


0001 IMPLICIT REAL*8 (A-H,O-Z) 

0002 DIMENSION DERY(8),Y(8),PRMT(5) ,AUX(8,8) 

0003 EXTERNAL FCT, OUTP 

0004 COMMON RFD, XD, XE, XAFD,RE,R,XFFD, XFFQ,XQ,RFQ, XFKD,XFKQ, 
1XKKD, XKKQ,V, XAKD, XAKQ, RKD, RKQ 

0005 COMMON WG leet eGo ice? Ol Se GoleC32 <3 3 ¢Dide Dl, Ol, Delle Dec. 
gD VAEICG DS Weed De Ve es en a tee, ab sp CARS BDV CRA Davy 7A 

0006 COMMON A(6,6),C(3,4) 

0007 COMMON 77, SEOLD, .REOLD, ZLOLD,VOLD, I; 11,J 

0008 COMMON /AA/VFD,VFQ,ET,BN,BK,E,S,Z3 

0009 BAD Cas dcl jis WGo CO) 144) 5 ¥ CO) VCO) 5 VY 0705. (8) 

0010 33 FORMAT(1P2D20.10/1P2D20.10/1P2D20.10/1P2D20.10) 

0011 READ (5, 30) RFD, XD, XE, XAFD, XAFQ, RE,R, XFFD, XFFQ,XQ, RFQ, XFKD, XFKQ, 
1XKKD, XKKQ,V, DD, VFD, VFQ, XAKD, XAKQ, RKD, RKQ 

COEZ P30 FORMAT CCsD10 5377000. 3/5010) 9/4P2D19).,.10/4D10. 3) 

0013 D0..2.K=1,,.8 

0014 2 DERY (K)=1.0d0/8.0D0 

0015 PRMT(1)=0. 

0016 PRMT (2)=1.5D0 

0017 PRMT (3)=0.001D0 

0018 PRMT (4)=0.001D0 

0019 NDIM=8 

0020 G=314.0D0/6.0D0*0.004D0 

0021 ZL=0.15D0/314.0D0 

0022 X1=XAFD 

0023 X2=XAKD 

0024 X3=-XAFQ 

0025 X4=-XAKQ 

0026 B13=R+RE 

0027 B14=XAFQ 

0028 B15=XAKQ 

0029 B16=- (XQ+XE) 

00 30 B21=-RFD 

0031 B32=-RKD 

00 32 F41=-XAFD 

0033 F42=-XAKD 

0034 F43=XD+XE 

0035 F46=R+RE 

0036 F54=-RFQ 

0037 F65=-RKQ 

0037 Al1=1.0D0/314.0D0*XAFD 

0039 Al2=1.0d0/314.0DO*XAKD 

0040 Al 3=-1.0D0/314.0DO0*XD 

0041 A21=1.0D0/314.0D0*XFFD 

0042 A22=1.0D0/ 314.0DO0*XFKD 

0043 A23=-1.0D0/314.0DO*XAFD 

0044 A31=1.0D0/314.0D0*XFKD 

0045 A32=1.0D0/ 314.0D0*XKKD 


0046 


A33=-1.0D0/314.O0D0*XAKD 
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AM=A11*A22%A33+A1 2 *A2 3%A31+A1 3%A21*A32-A31 *A22*A1 3- 


1A32%*A23%*A11-A33%*A21%*A12 
C11=(A22*A33-A32%*A23) /AM 
C12=(A32%*A1 3-A12*A33)/AM 
C13=(A12*A23-A22*A1 3) /AM 
C21=(A31*A2 3-A21*A33) /AM 
C22=(A11%*A33—A31*A13) /AM 
C2 3=(A21%*Al1 3-Al1*A2 3) /AM 
C31=(A21*A32-A31%*A22) /AM 
C32=(A31*A12-A11%*A32) /AM 
C33=(A11%*A22-A21%*A12) /AM 
E11=1.0D0/314.0D0*XAFQ 
E12=1.0D0/314.0D0*XAKQ 
E13=-1.0D0/314.0D0*XQ 
E21=1.0D0/314.0DO*XFFQ 
E22=1.0D0/314.0D0*XFKQ 
E23=-1.0D0/314.0D0*XAFQ 
E31=1.0D0/314.0D0*XFKQ 
E32=1.0D0/314.0D0*XKKQ 
E33=-1.0D0/414.0D0*XAKQ 


EM=E11*E22*E33+E12 *E2 3*E31+E1 3*E21*E32- 
1LE31*E22*E1 3—E32*E2 3%E11~-E33*E21*E12 


D11=(E22*E33-E32*E23) /EM 
D12=(E32*E1 3-E12*E33) /EM 
D1 3=(E12*E23-E22*E13) /EM 
D21=(E31*E23-E21*E33) /EM 
D22= (E11*E33-E31*E13) /EM 
D2 3=(E21*E13-E11*E2 3) /EM 
D31=(E21*E32-E31*E22) /EM 
D32= (E31 *E12-E11%*E32) /EM 
D33=(E1L1*E22-E21*E12) /EM 
A(1,1)=C11*B21 
A(1,2)=C13*B32 
A(1,3)=C11*B13 
A(1,4)=C11*B14 
A(1,5)=C11*B15 
A(1,6)=C11*B16 
A(2,1)=C22*B21 
A(2,2)=C23*B32 
A(2,3)=C21*B13 
A(2,4)=C21*B14 
A(2,5)=C21*B15 
A(2,6)=C21*B16 
A(3,1)=C32*B21 
A(3,2)=C33*B32 
A(3,3)=C31*B13 
A(3,4)=C31*B14 
A(3,5)=C31*B15 
A(3,6)=C31*B16 
A(4,1)=D11*F41 
A(4,2)=D11*F42 
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0096 A(4,3)=D11*F43 

0097 A(4,4)=D12*F54 

0098 A(4,5)=D13*F65 

0099 —A(4,6)=D11*F46 

0100 A(5,1)=D21*F41 

0101 A(5,2)=D21*F42 

0102 A(5,3)=D21*F43 

0103 A(5,4)=D22*F54 

0104 A(5,5)=D23*F65 

0105 A(5,6)=D21*F46 

0106 A(6,1)=D31*F41 

0107 A(6,2)=D31*F42 

0108 A(6, 3)=D31*F43 

0109 A(6,4)=D32*F54 

0110 A(6,5)=D33*F65 

0111 A(6,6)=D31*F46 

0112 1=0 

0113 CALL DRKGS (PRMT,Y,DERY,NDIM, IHLF,FCT,OUTP, AUX) 

0114 IF (IHLF.GT.10) WRITE (6,66) IHLF 

0115 66 FORMAT (1X,I3) 

0116 STOP 

0117 END 

0001 SUBROUTINE FCT(X, Y, DERY) 

0002 IMPLICIT REAL*8 (A-H,0-Z) 

0003 DIMENSION DERY(8) ,Y(8) ,PRMT(5) , AUX(8, 8) 

0004 COMMON RFD,XD, XE, XAFD, XAFQ,RE,R, XFFD, XFFQ,XQ,RFQ,XFKD, XFKQ, 
1XKKD , XKKQ,V, XAKD, XAKQ, RKD, RKQ 

0005 COMMONVCIT4'C12 C134 621, C22%0238, 6814 C825 COSh D114 D127 D1S),021 ,D22, 
ID23, Dole 304 Hl X24 RON RALDD, Gh Z1NZ23 18 z 

0006 A(6,6),C(3,4) 

0007 COMMON ZL, XEOLD,REOLD, ZLOLD, VOLD,I,I1I,J 

0008 COMMON /AA/VFD,VFQ,ET,BN,E,S,Z3 

0009 ED=RE*Y (3) -XE*Y (6) +ZL*DERY (3) +DSIN(Y8) ) *V 

0010 EQ=RE*Y (6)+XE*Y (3)+ZL*DERY (6)+DCOS(Y8) ) *V 

0011 ET= (ED¥*2+EQ**2) **0.5 

0012 IF(X.LE.0.001D0) GO TO 100 

0013 E=X1*Y (6) *A(1, 1) *Y¥(1)+X1¥*Y (6) ®A(1, 2) ¥Y (2) +X1*Y (6) *A(1, 3) *Y (3) + 


LX1*Y (6) *A(1, 4) *Y (4) +X1*Y (6) XA(1, SO*Y (5) +X1*Y (6) *A(1, 6) *Y (6) + 
LX1*Y (1) ¥A(6, 1) *¥(1)+X1*Y (1) *A(6, 2) *Y¥ (2) +X1*Y (1) *A(6, 3) *Y (3) + 
1XL*Y (1) *A (6,4) *Y¥(4)+X1*Y (1) *A(6, 5) *¥ (5) +X14*Y (1) *A(6, 6) *Y (6) + 
1X2*Y (6) ¥A(2, 1) *Y (1) +X2*Y (6) *A(2, 2) *Y¥ (2) +X2*Y (6) ¥A(2, 3) *Y (3) + 
LX2*Y (6) *A(2, 4) ¥Y (4) +X2*Y (6) *A(2, 5) *Y(5)+X24Y (6) ¥A(2,6) *Y (6) + 
1X2*Y (2) *A (6,1) *Y(1)+X2*Y (2) *A(6, 2) *Y (2) +X2*Y (2) ¥A(6, 3) *Y(3) + 
1X2*Y (2) *A (6,4) *Y (4) +X2*Y (2) *A(6, 5) *Y (5) +X2*Y (2) ¥A(6 6) *Y(6)+ 
1X3*Y (3) *A (4,1) *¥ (1) +X3*Y (3) *A(4, 2) *Y (2) +X3*Y (3) *A(4, 3) ¥Y(3) + 
1X3*Y (3) *A(4, 4) *Y (4) +X3*Y (3) *A(4, 5) *Y¥ (5) +X3*Y (3) ¥A(4 6) *Y (6) + 
1X3*Y (4) ¥A(3, 1) *Y (1) +X34Y (4) *A(3, 2) *Y (2) +X3*Y (4) ¥A(3, 3) *Y (3) + 
1X3*Y (4) *A(3, 4) *Y¥ (4) +X3*Y (4) *A(3, 5) *Y (5) +X3*Y (4) *A(3,6) *Y (6) + 
1X4*Y (3) ¥A(5,1) *Y¥(1)+X4*Y (3) *A(5, 2) *¥(2)+X4*Y (3) ¥A(5, 3) *Y (3) + 
1X4*Y (3) *A(5, 4) *Y (3) +X44*Y (3) *A(5,5) *Y (5) +X4*Y (3) *A(5 6) *Y (6) + 


} 
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LX4*Y (5) *A(3,1) *Y (1) +X4*Y (5) *A(3, 2) #Y (2) +X4*Y (5) *A(3, 3) *Y(3)+ 
LX4*Y (5) *A(3,4) *Y (4)+X4*Y (5) *A(3, 5) #Y(5)+X4*Y (5) *A(3,6) *Y (6)+ 
1 (XL*Y (6) *C11+X2 *Y (6) *C21+X3*Y (4) *C31+K4 *Y (5) *C31) *VADSIN(Y (8) )+ 
1 (XL *Y (1) #D31+X2*Y (2) *D31+X3*Y (3) *D11+X4*Y (3) *D21) *V*DCOS (Y (8) ) 
BK= (X1*Y (1) *D32+X2*Y (2) *D32+X3*Y (3) *D12+X4 *Y (3) *D22) *RFQ/XAFQ 
B1=DSIGN(1. ODO, BK) 

BN=(xL*Y (6) *C12+X2 *Y (6) *C22+X3*Y (4) *C32+X4*Y (5) *C32) *RFD/XAFD 
B2=DSIGN(1.0D0, BN) 

T1=-314.0D0/6.0D0* (E-5.0D0*B2*BN-5. ODO*B1*BK) 
T2=-314.0D0/6.0D0*(E+5. ODO*B2*BN+5 . OD0*B1*BK) 

Z1=Y (8)-1.0D0/ (G**2) *314. ODO*DERY (7) 

Z2=314.0D0*Y (7)+1.0D0/G*314.0DO*DERY (7) 

23=314. ODO*DERY (7) 

T=T1 

IF(Z2.GT.0.0) T=T2 

SS=-1. ODO/T*Z2*G**2 

IF(SS.GT.90.0D0) SS=90.0D0 

Z=G*Z3-T+T*DEXP (SS) 

IF(Z"GT?0.0) GO TO 71 

VFD=B2*5.0D0 

VFQ=B1*5.0D0 

GO TO 100 

VFD=-B2*5. ODO 

VFQ=-B1*5.0D0 


100 DERY(1)=A(1,1) *Y(1)+A(1, 2) *¥ (2) +A(1, 3) *¥ (3) +A(1, 4) *¥ (4) + 


1A(1,4) *Y(4) *Y (7) +A(1, 5) *¥(5)4+A(1, 5) *¥ (5) *¥ (7) +A(1, 6) *¥ (6) + 
2A(1,6) *Y(6) *¥ (7) +C11*V*DSIN(Y (8) )+C12*RFD/XAFD*VFD 

DERY (2) =A(2,1) *¥ (1)+A(2, 2) *¥ (2) +A(2, 3) *¥ (3) +A(2,4) *¥ (4) 4+ 
1A(2,4) *Y¥ (4) *¥ (7) +A (2,5) *¥ (5) *¥ (7) +A(2,6) *Y (6) +A(2,5) ¥Y(5)+ 
2A(2,6) *Y (6) *¥ (7) +C21*V*DSIN(Y (8) )+C22*RFD/ XAFD*VFD 

DERY (3)=A(3,1) *Y¥(1)+A(3, 2) *¥ (2) +A(3, 3) ¥¥ (3) +A(3, 4) *Y¥ (4) 4+ 
1A(3,4) *Y (4) *Y (7) +A(3, 5) ®¥ (5) +A(3,5) *¥ (5) ¥ (7) +A(3, 6) ¥Y (6) + 
2A(3,6) *Y (6) *¥ (7) +C31*V*DSIN(Y (8) )+C32*RFD/XAFD*VFD 

DERY (4) =A(4,1) *Y¥(1)+A(4, 1) *Y (1) *¥ (7) +A(4, 2) ¥¥ (2) + 

1A(4, 2) *¥ (2) ¥¥ (7) +A(4, 3) ®¥ (3) +A(4, 3) XY (3) ®Y¥ (7) +A(4, 4) *¥ (4) 4+ 
2A(4,5) *¥(5)+A(4, 6) *Y (6) +D11*V*DCOS (y (8) )+D12*RFQ/XAFQ*BFQ 
DERY (5) =A(5,1) *¥ (1) +A(5,1) *Y(1) *¥ (7) +A(5, 2) *¥ (2) + 

1A(5,2) *¥ (2) *¥ (7) +A(5, 3) *¥(3)+A(5, 3) *Y¥ (3) *¥ (7) +A(5,4) *¥ (4) + 
2A(5,5) *Y¥(5)+A(5, 6) *¥ (6)+D21*V*DCOS (Y (8) )+D22*RFQ*VFQ/XAFQ 
DERY (6) =A(6,1) *Y¥(1)+A(6,1) *Y(1) *¥ (7) +A(6, 2) *¥(2)+ 
1A(6,2) *Y¥(2) *Y¥(7)+A(6, 3) *¥ (3) +A(6, 3) *¥ (3) *¥ (7) +A(6, 4) *¥ (4) + 
2A(6,5) *Y¥(5)+A(6,6) *Y (6)+D31*V*DCOS (Y (8) )+D32*RFQ*VFQ/XAFQ 
IF(X.<LT.0.06D0) GO TO 10 

DERY (7)=-1.0D0/6.0DO* (XAFD*Y (1) *Y (6) +XAKD*Y (2) *Y¥ (6) - 
LXAFQ*Y (4) *Y (3) -XAKQ8Y (5) *¥(3) )- 
21.0D0/6.0D0*DD*314.0D0*Y (7)+ 
31.0D0/6.0D0*1.051999999999999*1.0D0 

CONTOS25 

DERY (7)=-1.0D0/6. ODO* (XAFD*Y (1) *Y (6) +XAKD*Y (2) *Y (6) - 
LXAFQ*Y (4) *Y¥ (3) -KAKQ*Y (5) *¥ (3) )- 
21.0D0/6.0DO*DD*314.O0DO*Y(7)+ 
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31.0D0/6.0D0*1.051999999999999D0*2. 0D0 
25 DERY(8)=314.0D0*Y (7) 

RETURN 

END 


SUBROUTINE OUTP (X,Y, DERY, IHLF ,NDIM, PRMT) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION DERY(8),Y(8) ,PRMT(5) ,AUX(8,8) 

COMMON RFD, XD, XE,XAFD,XAFQ,RE,R,XFFD,XFFQ,XQ,RFQ, XFKD, XFKQ, 
1XKKD, XKKQ,V, XAKD, XAKQ, RKD, RKQ, 

COMMON: C114 C12;C13°C21022 023,031, C32,033,D11,D1L2501 3.0216 p22, 
PN25eD 31S D324 Dos KPO ROE XC DDI Gs 74 p22 Tez 

COMMON A(6,6),C(3,4) 

COMMON ZL,XEOLD, REOLD,ZLOLD,VOLD,I,II,J 
COMMON/AA/VFD,VFQ,ET, BN, BK,E,S,Z3 

IF(I/4*4.EQ.1) WRITE(6, 88)X,VFD,VFQ,Z,Y(1),Y(3),Y(4) ,ET,Y(7), 
1DERY (7) ,Y (8) 

88 FORMAT(1X,1P10D11.3,1P1D13.5) 

I=I+1 

RETURN 

END 
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